The partition representation of enzymatic reaction networks and its application for searching bi-stable reaction systems

The signal transduction system, which is known as a regulatory mechanism for biochemical reaction systems in the cell, has been the subject of intensive research in recent years, and its design methods have become necessary from the viewpoint of synthetic biology. We proposed the partition representation of enzymatic reaction networks consisting of post-translational modification reactions such as phosphorylation, which is an important basic component of signal transduction systems, and attempted to find enzymatic reaction networks with bistability to demonstrate the effectiveness of the proposed representation method. The partition modifiers can be naturally introduced into the partition representation of enzymatic reaction networks when applied to search. By randomly applying the partition modifiers as appropriate, we searched for bistable and resettable enzymatic reaction networks consisting of four post-translational modification reactions. The proposed search algorithm worked well and we were able to find various bistable enzymatic reaction networks, including a typical bistable enzymatic reaction network with positive auto-feedbacks and mutually negative regulations. Since the search algorithm is divided into an evaluation function specific to the characteristics of the enzymatic reaction network to be searched and an independent algorithm part, it may be applied to search for dynamic properties such as biochemical adaptation, the ability to reset the biochemical state after responding to a stimulus, by replacing the evaluation function with one for other characteristics.


Introduction
The intracellular signal transduction system functions as a mechanism to control cell proliferation, apoptosis, differentiation, and homeostasis, and its malfunction is thought to be one of the substantial causes of cancer formation for cells. The regulatory mechanism of signal transduction systems is closely related to the mechanism of action of drugs, such as anticancer drugs, and has therefore been the subject of extensive research in recent years. In addition, from the viewpoint of synthetic biology, which has been the focus of much research in recent years, design methods for signal transduction systems have become necessary. However, because biochemical reactions are nonlinear, it is difficult to establish theoretical analysis methods, and signal transduction systems in particular have not been systematically analyzed because their interactions are more complex than those of intracellular metabolic systems, which they control. The current situation is that the parameters are fixed for each system to be analyzed, and the analysis is performed by computer simulation. The MAPK cascade, a representative signal transduction system, which relays cell growth factor (EGF) from the cell membrane to the cell nucleus, and its abnormality is thought to be the cause of cell canceration, and much knowledge has been obtained [1][2][3][4][5][6]. A major component of the signal transduction system represented by the MAPK cascade is the cyclic reaction system, which activates and deactivates enzymes by phosphorylation and dephosphorylation. The cyclic reaction system is a combination of two post-translational modification reactions. Therefore, we have been conducting a comprehensive analysis of the stability of regulatory networks that consist of activation and inactivation cyclic reaction systems of enzymes and their mutual controls [7]. There, the control relationship between the cyclic reaction systems is represented by a control matrix. The control matrix, which is an adjacency matrix, has 1 as its i-th row and j-th column components when the activating enzyme of cyclic reaction system j catalyzes the activation of cycle reaction system i, and -1 when it catalyzes the deactivation.
In a similar study, Kuwahara et al. [8] exhaustively analyzed the influence of the control structure on the stochastic properties for regulatory networks of three to five nodes with only one feedback control. They also mentioned the control structures in which bistability emerges. Ramakrishnan et al. [9], Shah et al. [10], and Siegal-Gaskins et al. [11] have conducted exhaustive analyses of the ultrasensitive properties and bistability of signaling systems. Ultrasensitivity was reported as a response characteristic of the MAPK cascade, a typical signaling system mentioned above [1]. In a system with ultrasensitivity, the steady-state enzyme concentration of the output chemical species with respect to the enzyme concentration of the input chemical species increases rapidly after a certain threshold of the input. Bistability is when this ultrasensitivity is more strongly expressed, and hysteresis appears in the change of output to input. Thus, in a region with input, two steady state values of the output appear. In particular, it is said to be resettable if it has the property of being able to mutually transition between two steady states by varying the total concentration of the input chemical species within the range of protein concentrations expected in the cell [12].
Ma et al. [13] and Yao et al. [14] analyzed the biochemical adaptation and robustness properties using a similar control network of three nodes. Biochemical adaptation is the property that when the concentration of an input chemical species is increased in a steady-state system, the response of the output chemical species transiently increases and then decreases to a value close to the original steady-state value [13]. Robustness is the property that the various response properties of a system are less sensitive to the concentration of other chemical species in the system [12].
Adler et al. [15] focused on fold-change detection, which is a dynamic property in which the temporal pattern of the system's output in response to a transient input change in the cell depends on the ratio rather than the difference of the input changes, demonstrating the effectiveness of comprehensive analysis using regulatory networks. In these studies, the enzymatic reaction mechanism is mainly based on the Michaelis-Menten approximation or a simplified linear equation in order to reduce the computational complexity. However, it has been reported that approximations can lead to inaccuracies in the bistability and dynamics of the system [16,17]. In order to accurately analyze the properties of a system, it would be better to describe the system using only the law of mass action, but this would increase the number of parameters of the system and make the originally huge search space even larger, making exhaustive analysis difficult to carry out. Therefore, there is a need to develop a method to search for enzymatic reaction networks with useful properties such as ultrasensitivity, bistability, biochemical adaptation, robustness, and fold-change detection as mentioned above [18][19][20].
In this study, we focus on the enzymatic reaction network, which is formulated as a regulatory network in which post-translational modification reactions are the elements and they mutually regulate each other. Post-translational modification reactions, such as phosphorylation, are smaller units than the cyclic reaction systems. Instead of the control matrix described above, we propose to use a partition of the set as its representation, which is considered to be more suitable for search. In addition, we report the results of our search for a bistable and resettable enzymatic reaction network with the aim of demonstrating its effectiveness. Resettable is the property of being able to transition between two steady states by changing the total concentration of the input chemical species within the range of expected protein concentrations in the cell [12].

Partition representations for enzymatic reaction networks
A post-translational modification reaction with Michaelis-Menten-type enzymatic reaction as the reaction mechanism is considered as the basic component, considering to construct the enzymatic reaction network composed of N components. In this regard, reactions in which multiple enzymes bind to a single substrate simultaneously are not considered for simplicity in this formulation. The Michaelis-Menten-type enzymatic reaction is a series of chemical reactions in which an enzyme (E) binds to a substrate (S) to form a temporary enzyme-substrate complex (C), the substrate is transformed into a product (P) on the complex, and then decomposes into the product (P) and the enzyme (E), as shown in Eq (1), where the reaction rate constants are denoted by a, d, and k.
The substrates, enzymes, and products of N Michaelis-Menten-type enzymatic reactions are denoted by S i , E i , and P i , respectively, and the set is denoted by M. That is, M = {S 1 , E 1 , P 1 , S 2 , E 2 , P 2 , . . ., S N , E N , P N }. Take one partition P of the set M. A partition is a family of mutually disjoint subsets of a set, and the union of those subsets is the original set. By identifying the elements of the partition, i.e., the chemical species that belong to one subset of the original set M, an enzymatic reaction network is constructed from multiple post-translational modification reactions. Fig 1 shows an example of a basic enzymatic reaction network consisting of post-translational modification reactions. The primary building blocks of post-translational modification reactions are numbered consecutively; the substrate, enzyme, and product of the i-th post-translational modification reaction are labeled S i , E i , and P i , respectively. If the chemical species to be equated is the substrate S i and the product P i , the names of the species to be equated are listed in the rectangle. The enzyme E i is represented by the small circle between the red arrows, and if there is a chemical species that is identical to the enzyme, it is connected to that chemical species by a dotted line. It may be easier to understand intuitively if you think of the small circle as representing the enzyme-substrate complex. The identification of P i and S j is a representation of the two-step enzymatic reaction from S i to P j . Fig 1C and 1D 1 , P 2 }} of M, respectively. The identification of S i and S j , and P i and P j , respectively, is a representation of the branching into and confluence from two enzymatic reactions.  Examples of enzymatic reaction networks composed of posttranslational modification reactions as the primary building blocks; a: product of enzymatic reaction 1 catalyzes enzymatic reaction 2, b: two-step reaction, c: branching reaction, d: confluence reaction, e: autocatalysis, f: association-dissociation reaction, g: dimer formation separation reaction. The primary building blocks of posttranslational modification reactions are numbered consecutively; the substrate, enzyme, and product of the i-th post-translational modification reaction are labeled S i , E i , and P i , respectively. If the chemical species to be equated is the substrate S i and the product P i , the names of the species to be equated are listed in the rectangle. The enzyme E i is represented by the small circle between the red arrows, and if there is a chemical species that identifies with the enzyme, it is connected to that species by a dotted line. and product and the enzyme S 1 = P 1 = E 1 associate and dissociate themselves, corresponding to the dimer formation separation reaction. The identification of S i, P i , and E i is an expression of the dimer formation-separation reaction. Fig 2 shows examples of more complex enzymatic reaction networks and their partition representations. Fig 2A shows an example of an enzymatic reaction network consisting of N = 4, i.e., four post-translational modification reactions, where the original set of the partition is M = {S 1 , E 1 , P 1 , S 2 , E 2 , P 2 , S 3 , E 3 , P 3 , S 4 , E 4 , P 4 } and the partitions is {{S 1 , P 2 }, {P 1 , S 2 , E 1 , E 4 }, {S 3 , P 4 }, {P 3 , S 4 , E 2 , E 3 }}. This is an example of an enzymatic reaction network consisting of two cyclic reaction systems with auto-activating and mutually inhibitory feedbacks, which is an typical of enzymatic reaction network with bistable properties. The MAPK cascade, one of the representative signaling systems, is shown in Fig 2B. The MAPK cascade consists of a cascade of three kinases, MAPKKK, MAPKKK, and MAPK. The activated enzyme activates the phosphotransferase in the next step, thereby transmitting signals within the cell. An example of an enzymatic reaction network consisting of 10 post-translational modification reactions, with M = {S 1 , E 1 , P 1 , S 2 , E 2 , P 2 , . . ., S 10 , E 10 , P 10 10 , S 7 }, {P 7 , P 9 , S 8 , S 10 }, and {P 8 , S 9 } are inactive, one-phosphorylated, and two-phosphorylated MAPKs, respectively.
The advantage of the partition representation is that every partition corresponds to an enzymatic reaction network, and thus there is a one-to-one relationship between the partition and the enzymatic reaction network. The number of different enzymatic reaction networks consisting of N post-translational modification reactions is the same as the total number of different partitions of the set with number of elements 3N, a number whose value is known as the Bell number B 3N , expressed in the recurrence equation as in Eq (2), where B 0 = B 1 = 1.
For example, the total number of enzymatic reaction networks consisting of four post-translational modification reactions, as shown in Fig 2A, The symbol is a function that returns the number of elements in a set. Therefore, is the number of elements that are not common to the sets a i and b j . For example, let the source set of the partition be It can be seen that the similarity of the enzyme reaction network correlates with this distance. Furthermore, this distance between partitions satisfies the triangle inequality and can be used to visualize the distribution of partitions.

Derivation of differential equation systems and conservation laws from partitions
From the partition, we can derive a system of differential equations and conservation laws that describe the behavior of the corresponding enzymatic reaction network as follows. If the number of post-translational modification reactions constituting the enzymatic reaction network is N, the total number of elements in the set from which the partitioning is made is 3N, since three elements, that is, a substrate, an enzyme, and a product, correspond to one post-translational modification reaction. Therefore, if we denote the elements by x i , the source set of the partition can be expressed as Let one of the partitions of M be A = {a 1 , a 2 , . . ., a m }, and let s k be the set of subscript i of x i which is an element of the kth element a k of the partitions A. For this partition A, a system of differential equations can be generated by the following procedure.
1. Generate 4N reaction rate equations for substrate S i , enzyme E i , enzyme-substrate complex C i , and product P i in each post-translational modification reaction, as shown in Eq (4).
2. For each element a k in the partition, generate a chemical species y k that identifies the all elements of a k , and replace the original reaction rate equations for the elements of a k with the following Eq (5), while all elements of a k appearing on the right side of Eq (5) are replaced by y k . The s k is the subscript set of the elements of a k .
For the derivation of Eq (4), if only the law of mass action is used as the reaction equation for the post-translational modification reaction, and the variable name is expressed as x[chemical species], the reaction equation derived is the following for the post-translational modification reaction i.
For example, if the original set of partitions is M = {S 1 , E 1 , P 1 , S 2 , E 2 , P 2 }, and the partition is A = {{S 1 }, {E 1 }, {P 1 , E 2 }, {S 2 }, {P 2 }} as seen in Fig 1A, then by applying the law of mass action, the above Eq (4) becomes the following Eq (7), where the variable name is represented by x The equations corresponding to Eq (5), obtained by equating P 1 and E 2 according to division A, are Eq (8). In particular, when some chemical species are considered identical, they are listed in brackets.
For the conservation law, we first generate a set of lists Q = {q 1 , . . ., q i , . . ., q 2N } of chemical species whose total concentration is conserved for each post-translational modification reaction. The reason why the total number is 2N is that for each post-translational modification reaction, there is a conservation law for the enzyme and a conservation law for the substrate. Next, the following procedure is applied sequentially to each element a k of the partition A, in turn.
1. Concatenate all the list of elements in Q that contain elements in common with a k , and let r k be the list.
2. Remove the elements of a k from the list rk and add the chemical species y k , which is identical to the elements of a k .
3. Add r k to the element of Q that does not contain an element in common with a k , and make it Q again.
It is important to note that the above procedure concatenates, rather than sums, the elements of Q that contain elements in common with a k . This is because the overlap is intrinsic when the list to be combined contains the common enzyme-substrate complex C i .
If only the law of mass action is used as the reaction equation for the post-translational modification reaction, and the variable name is represented by x[chemical species], the initial elements of the set of conservation laws Q to be derived are the following for the post-translational modification reaction i.
For example, for the enzymatic reaction networks in Fig 1A,

Resettabe bistability
A search for bistable and resettable enzymatic reaction networks was attempted to demonstrate the effectiveness of the partition representation of enzymatic reaction networks. Resettable is the property of being able to mutually transition between two steady states by varying the total concentration of the input chemical species within the range of expected protein concentrations in the cell [12]. A typical resettable bistable enzymatic reaction network is shown in Fig 2A. The relationship between bistability and resettability when the input chemical species is E 1 and the output chemical species is P 2 is shown in Fig 4. For example, in the enzymatic reaction network of Fig 2A, the input chemical species E 1 corresponds to the small circle labeled 1, and the output chemical species P 2 corresponds to the rectangle in the lower left corner that equates S 1 and P 2 . The horizontal axis is the logarithm of the input species concentration with a base of 2, and the vertical axis is the total output species concentration normalized to between 0 and 1. The range of protein concentrations expected in the cell is 2 −5 to 2 5 in m mol/m 3 .
The curve represents the concentration of the output species at equilibrium, where the solid line corresponds to the stable equilibrium point, or steady state value, and the dashed line corresponds to the unstable equilibrium point. Fig 4A shows the resettable bistable aspect. The solid blue curves at both ends of the graph correspond to a single steady state, which are monostable states. The red curve in the center has two steady states, corresponding to bistable states. As the total concentration of the input species is increased from the smaller one, the black arrow on the right makes a discontinuous jump from the upper steady state to the blue monostable state on the right. When the total input species concentration is reduced from that state, the system jumps to the steady state below at the black arrow on the left side. Thus, it is a resettable bistable because it can mutually transition between two steady states. Fig 4B is an example where the curve in Fig 4A extends to the right and the region of monostability on the right exceeds the range of protein concentrations expected in the cell. In this case, the black arrows show that it is possible to jump from the upper steady state to the lower steady state, but it is not possible to jump from the lower steady state to the upper steady state, even if the total input chemical species concentration is increased within the range of expected protein concentrations in the cell. This is an example of a bistable that is not resettable. Fig 4C shows an example where the graph extends further to the left. In such a case, if the initial state of the system is above the red dotted line in the center, it shifts to the upper steady state, and if it is below, it shifts to the lower steady state, and the state cannot be changed even if the total concentration of the input chemical species is changed within the range of the expected protein concentration in the cell. This is another example of bistability that is not resettable.

Exploration of enzymatic reaction networks in the partition representation space
For the partition representation of the enzymatic reaction network, an algorithm to search for the enzymatic reaction network with the desired response characteristics can be implemented in a natural way by sequentially applying the following two partition modifiers. • Join partition modifier: Combines two randomly selected elements. For example, in the example of Fig 3 above, by combining the second and third elements of the partition corresponding to the enzymatic reaction network A, the partition corresponding to the enzymatic reaction network C is obtained.
First, we set up an evaluation function ϕ(P,ν) that quantifies the desired response properties of the dynamics or steady state of the enzymatic reaction network generated from the partition. ϕ(P,ν) is a function of the partition P, with the parameter ν, which is a pair of reaction rate constants of each post-translational modification reaction and the total concentration of enzymes that constitute the enzymatic reaction network corresponding to the partition.
In the following algorithm, we use the function Ф(P) to perform a random search on the parameter value ν for a given partition P. The function Ф(P) quantifies the evaluation value at λ randomly chosen parameter values from the set of candidate values Λ of the parameter ν by ϕ(P,ν) and returns the pair of the largest evaluation value and the parameter ν at that time as the value.

PLOS ONE
The partition representation of enzymatic reaction networks The main loop of the algorithm is as follows.
1. Determine one partition P at random.
2. If the value of the evaluation function Ф(P) is less than ρ or the number of iterations is less than γ, repeat P = ψ(P, P, σ).
3. Return P as the value.
The partition search function ψ(P 0 , P, σ) is a recursive function whose arguments are the initial partition P 0 , the partition P, and the partition depth σ described below. ψ(P 0 , P, σ) returns the partition whose evaluation function is greater than or equal to Ф(P 0 ).
1. If σ is zero, return P 0 as the value.
2. Generate P' from P with the join partition modifier.
10. If none of the above conditions are satisfied, return P 0 as the value.
We apply this random search algorithm to the aforementioned search for enzymatic reaction networks with resettable bistability with the aim of demonstrating the effectiveness of the partition representation of enzymaric reaction networks. In other words, we design an evaluation function ϕ(P,ν) for resettable bistability and use it in this algorithm.

Evaluation function for resettable bistability
Bistability is the property of having two steady state values for a given parameter value ν, as described above. Furthermore, resettable is the property that allows the transition between these two steady states by changing the value of the total concentration of the enzyme containing the input chemical species. Here, we create an evaluation function ϕ(P,ν) that gives a high evaluation value for a typical resettable bistable aspect as shown in Fig 4A, i.e., an enzymatic reaction network that is monostable at the edge of the range taken by the total concentration of input chemical species and bistable in the central part.
ϕ(P,ν) is a function of the partition P, with the reaction rate constants of each post-translational modification reaction and the total concentration of enzymes comprising the enzymatic reaction network corresponding to the partition as parameters ν. The variations in the total concentration of enzymes containing the input chemical species were set to 21 discrete values of O = 2 −5 , 2 −4.5 , . . ., 2 5 . These values, as well as the reaction rate constants to be set as ν, are set to include approximately the values reported for the kinases that make up the MAPK cascade, which is known to be a typical cellular signaling system.
The following procedure corresponds to the aforementioned evaluation function ϕ(P,ν). P and ν are given as arguments.
1. For a given parameter value ν, only the total concentration of the enzyme containing the input chemical species is varied sequentially and discretely for 21 values of O = 2 −5 , 2 −4.5 , . . ., 2 5 , the total concentration of each enzyme is randomly assigned as the initial value so as to satisfy the conservation law of each enzyme. The steady state values were obtained by numerically solving the generated differential equations and making convergence judgments as appropriate. The details are described in the next section. Then 21 pairs of steady state values of the output chemical species are obtained. Here, the steady-state value of the output enzyme species is normalized to its total concentration, and is a value between 0 and 1.

This is repeated
Resettable bistability is the property of being monostable at the edge of the range taken by the total concentration of the input chemical species and bistable in the middle part, as explained in Fig 4. The eVL evaluates the degree of this property. Its input V is a vector of the above 21 quadruple variance values and V i denotes the each element. The value of the function eVL was normalized by a maximum value so that it would be between 0 and 1. The α is the coefficient for this purpose, which is the reciprocal of the maximum value 0.604762 that the summation equation on the right side takes when V = {0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0}. The 3D shape of the function wVL that makes up the function eVL is shown in Fig 5. Its input v is the quadruple variance value V i and l is the position i of V i in the vector V. The function wVL is designed to be maximal for perfect monostability (v = 0) at the edge of the region, and also to be maximal for perfect bistability (v = 1) at the center of the region, and continuous intermediate evaluation values in the middle of the region.

Hyperparameters of the search algorithm
The following is a summary of the hyperparameters included in the algorithm for searching for resettable bistable enzymatic reaction networks, which is an application of the partition representation, and the values used in the actual search. First, the hyperparameters that are independent of the characteristics of the enzymatic reaction networks to be explored include the following. As mentioned earlier, the enzymatic reaction network shown in Fig 2B is known as a typical resettable bistable enzymatic reaction network [7], and the number of post-translational modification reactions that make up the network is 4, so the value of N was set to 4. The input and output chemical species are set to be E 1 and P 2 , respectively, which belongs to the different post-translational modification reactions, because the enzymatic reaction network explored would be severely limited if the input and output chemical species were those of the same posttranslational modification reaction.
We did not know whether the values we chose for the depth of recursive search σ and the maximum number of loops γ were optimal or not. After trying several combinations, we chose the one that would allow us to explore a single enzyme reaction network in at least a few days. The set of values Λ used for the total concentration of enzymes and the rate constants of the enzymatic reactions was set to include approximately the values in mmol/m 3 -s system reported for the phosphatases comprising the MAPK cascade [21][22][23][24][25]. For λ, we tried several other values, but did not find much difference. Hyperparameters specific to the search target, i.e., the search for bistable enzyme reaction networks, include the following. In general, these will alter for different search targets.
• A series of total concentrations of the input chemical species O: {2 −5 , 2 −4.5 , . . ., 2 5 } • Number of steady-state values with random initial values θ: 5 We used the same variation of Λ as described above for O. We tried several other values, but did not find much difference for θ.
To obtain the above steady state values, the generated differential equations were solved numerically many times while updating the initial values, and the convergence was judged. This part has the following hyperparameters, which are common when steady state values are needed in the search process. Wolfram Mathematica's NDSolve function was used to solve the differential equations numerically. The protein concentration after τ seconds was calculated by NDSolve, and was judged to be converged when all of the relative error of the change, i.e., the ratio of the change to the total concentration of the protein, was less than ε. If there is no convergence, the protein concentration after another τ seconds is calculated and the same decision is made. If it does not converge after repeating this process κ times, it returns the result that it did not converge. The product of τ and κ, 3600 seconds, or 1 hour, was taken as the critical time to reach steady state. This time was taken as a reference for the response time of the signaling system in the cell. The values of the calculation time τ for one step and the convergence decision error ε were adopted after trial and error.

Enzymatic reaction networks found by the search algorithm
The search was conducted on a machine with a 3.4GHz i7-3770 processor CPU, 16GB of memory, and a Windows 10 operating system. The search program was implemented in Wolfram Mathematica 12.0. The bistable enzymatic reaction networks found in the search at the values of the hyperparameters described above are shown in Fig 6A-6H. Each of them uses a different seed of random numbers for the initial setting. The maximum CPU time required to discover a single enzymatic reaction network was approximately 30 hours. For example, the enzymatic reaction network in Fig 6A consists of three chemical species: the input is the chemical species that identifies E 1 , E 2 , P 1 , P 4 , S 1 , and S 4 , and the output is the chemical species that identifies P 2 and S 3 . The rest is a chemical species that identify E 3 , E 4 , P 3 , and S 2 . In particular, the input chemical species is in equilibrium with its dimer. Fig 6G is identical to the typical bistable enzymatic reaction network shown in Fig 2A. The reaction rate constants for each enzymatic reaction network, the total concentration of each chemical species, and the derived differential equation systems are shown in S1-S4 Tables provided as the supporting information, respectively. Fig 7 shows the bistable aspect of each enzymatic reaction network, where a-h corresponds to a-h in Fig 6. In order to plot the bistable aspect, the Monte Carlo method was used to find the steady state values of the output chemical species for each value of the total concentration of the enzyme containing the input chemical species. That is, when numerically solving the system of differential equations derived from the partition, the initial values were randomly allocated to the chemical species included in each conservation law. We can see that all the  obtained enzymatic reaction networks exhibit resettable bistability. Resettable is the property of being able to go back and forth between two bistable steady state values by changing the value of the input.   Fig 8A shows the evolution of the evaluation value as the search progresses. The horizontal axis is the number of times the structure or parameter values of the enzymatic reaction network were changed. The vertical axis is the evaluation value, with the maximum being 1. The legends a-h corresponds to a-h in Fig 6. Since the evaluation value reaches 0.6 in all cases in the first few steps, and almost all graphs overlap up to that point, the range of the vertical axis is set to 0.55 or higher so that there is less overlap thereafter. The search is limited to 150 times, and the search is terminated when the evaluation value exceeds 0.8. It can be seen that the aspect of convergence varies depending on the seeds of random numbers used at initial settings. Fig 6B shows movements within the partition space,

PLOS ONE
The partition representation of enzymatic reaction networks based on the distance between the partitions defined by Eq (3). Wolfram Mathematica's Graph function was used for drawing. The distance between all partitions was measured by the introduced distance function and the values were used in the EdgeWeight option. The GrpahLayout option specifies SpringEmbedding, which is arranged so that the total energy of the mutually coupled partitions is minimized by a spring of force equivalent to the distance between the partitions. The legends a-h corresponds to a-h in Fig 6. The red circles denote the final partitions in each search. We can see that the search is moving through the entire partition space.

Discussion
We have shown that an enzymatic reaction network composed of post-translational modification reactions can be represented by a partition of the set whose elements are the enzymes, substrates, and products. By identifying elements of each subset in the partition as the same chemical species, an enzymatic reaction network can be constructed. In particular, by equating the substrate and product within the same post-translational modification reaction, we can also describe the association and dissociation reaction between enzyme and substrate (product). However, it does not represent the binding of multiple proteins as found in scaffold proteins and as seen in membrane receptor proteins. As an extension of the partition representation, enzyme-substrate complexes may be added as elements of the original set, but an algorithm for deriving the system of differential equations and conservation laws from the partition needs to be devised.
Furthermore, the partition representation of enzymatic reaction networks was applied to the search for bistable networks. The two partition modifiers introduced into the partition representation worked well and various bistable enzymatic reaction networks were obtained. One of the enzymatic reaction networks obtained in the search were the typical bistable enzymatic reaction network in which two cyclic reaction systems with positive auto-regulatory feedbacks mutually negatively regulate each other.
Regarding the performance of the search algorithm, it was confirmed that the search was completed within an acceptable processing time. It is especially noteworthy that the processing time required for the search was greatly reduced compared to an exhaustive search. In addition, the visualization of the search path using the distance between the partitions introduced in this study was able to show the aspect of the search.
It is not difficult to increase the granularity of the primary building blocks from post-translational modification reactions to cyclic reaction systems. Instead of the original set of elements being the substrate, enzyme, and product of the post-translational modification reaction, we can have four types of enzymes: the active enzyme, the inactive enzyme, the activating enzyme that catalyzes the reaction that turns the inactive enzyme into the active enzyme, and vice versa. We plan to try it in the future.
For the hyperparameter of the search algorithm, which is the search depth, and the parameter values for evaluating the structure of an enzymatic reaction network, which are the reaction rate constant and the number of times the total concentration of each chemical species is randomly tested, we set the search depth to 4 and the number of random tests to 5, but we did not systematically test other values. It is possible that adjusting these values can speed up the process. Further speed-up can be achieved by using GPUs or by using tQCCM [16], an improved form of the Michaelis-Menten approximation, as a formulation of the reaction mechanism instead of the law of mass action.
Although the bistability targeted in this study is a steady-state property, the search algorithm is divided into an evaluation function specific to the characteristics of the enzymatic reaction network to be explored and an independent algorithmic part, so that it may be possible to apply it to the exploration of dynamic properties such as biochemical adaptation, by replacing the evaluation function. We plan to try it in the future. The number of post-translational modification reactions, which is the primary building block, was set to be 4, but even when the number was set to be 3, a bistable enzymatic reaction networks could be found. If the mechanism for automatic increase or decrease in the number of post-translational modification reactions can be included in the search algorithm, it will be possible to find a more optimal enzyme reaction network.  Table. The total concentration of enzymes in each enzymatic reaction network found. Each line corresponds to a-h in Fig 6. The second column is a list of sets of chemical species corresponding to the conservation laws. The third column is their total concentration in the same order as the second column. However, for the set containing E 1 , the input chemical species, we have moved it in the range of 2 -5 to 2 5 , instead of the values in the table. (XLSX) S3 Table. Differential equation system for each enzymatic reaction network found (a-d).

Supporting information
The first column corresponds to a-d in Fig 6. The second column is the name of the variable to be differentiated on the left side of the differential equation and the third column is the right side of the differential equation. The variable name is represented by x[chemical species name]. In particular, when some chemical species are considered identical, they are listed in brackets. (XLSX)

S4 Table. Differential equation system for each enzymatic reaction network found (e-h).
The first column corresponds to e-h in Fig 6. The second column is the name of the variable to be differentiated on the left side of the differential equation and the third column is the right side of the differential equation. The variable name is represented by x[chemical species name]. In particular, when some chemical species are considered identical, they are listed in brackets. (XLSX) S1 Data. (ZIP)